%tsData - GNSS station time series data object
% 
%   tsDataOut = tsData
%   tsDataOut = tsData(FILENA)
%
% The output tsData object is generated with the data taken from a .tenv3, 
% .tenv, .kenv or .pos ASCII file FILENA (the function that extracts these
% data is a method of this object). If FILENA is undefined or empty, the 
% filename can be interactively managed by means of a combo box.
%
% The properties of the generated object of tsData class are related to:
%   1) input time series (in this case, the data are managed by means of
%      tsData methods);
%   2) time series processing (in this case, no tsData methods act on the 
%      data).
%
% Complete list of properties 1) (please note that all these time series
% properties are defined only in the case of data taken from a tenv3 file. 
% In the case of tenv or kenv files some properties are undefined and, 
% therefore, the corresponding values are empty): 
%
%   Scalars/strings:
%
%         statName    : station name (string)
%         statLat     : station latitude (single value, degrees)
%         statLon     : station longitude (single value, degrees)
%         statHeight  : station heigh (single value, m)
%
%   Arrays:
%
%         dateS       : date (string YYMMMDD) 
%         dateyfrac   : date (fractional year)
%         MJD         : modified julian date  
%         GPSweek     : GPS week
%         GPSday      : GPS day
%         t           : date (MATLAB serial form)
%         reflon      : reference meridian longitude (degrees)
%         E0          : easting (m), integer portion (from ref. meridian)
%         Ed          : easting (m), fractional portion
%         E           : easting (m), complete
%         N0          : northing (m), integer portion (from equator)    
%         Nd          : northing (m), fractional portion
%         N           : northing (m), complete
%         V0          : vertical (m), integer portion (from equator)    
%         Vd          : vertical (m), fractional portion
%         V           : vertical (m), complete
%         antH        : antenna height (m) assumed from Rinex header
%         sE          : east sigma (m)
%         sN          : north sigma (m)
%         sV          : vertical sigma (m)
%         cEN         : east-north correlation coefficient
%         cEV         : east-vertical correlation coefficient
%         cNV         : north-vertical correlation coefficient
%
% Complete list of properties 2), which are related to the time series 
% processing and, therefore, are empty before such a processing:
%
%         cleanedE      : cleaned time series after outlier removal - East
%         cleanedN      : cleaned time series after outlier removal - North
%         cleanedV      : cleaned time series after outlier removal - Vertical
%
%         offsetsE      : offsets East
%         offsetsN      : offsets North
%         offsetsV      : offsets vertical
%         offsetsOpts   : offsets search general options
%   The value of each property offsetsE, offsetsN and offsetsV, as the 
%   offset search is carried out, is an 'offsets' struct variable (see 
%   below). The value of offsetsOpts:
%       - if the offsets recognition is carried out in an automatic way,
%         it is OptsGen.Offsets, where OptsGen is the struct variable with 
%         the options for StaVelMain computations;
%       - if the offsets recognition is carried out in the manual way,
%         it is the string 'Manual offsets recognition';
%       - if the offsets are taken from NGL database, it is the string 
%         'Offsets taken from Nevada Geodetic Laboratory database';
%       - if the offsets are taken from another database, it is the string
%         'Offsets taken from file ....' 
%
%         outliersE     : outliers East
%         outliersN     : outliers North 
%         outliersV     : outliers Vertical
%         outliersOpts  : outliers search general options
%   The value of each property outliersE, outliersN, outliersV, as the 
%   outlier search is carried out, is an 'outliers' struct variable (see 
%   below). The value of outliersOpts is OptsGen.Outliers, where OptsGen
%   is as above.
%
%     estimatedTrendE   : estimated trend East 
%     estimatedTrendN   : estimated trend North 
%     estimatedTrendV   : estimated trend Vertical
%     estimatedTrendOpts: estimated trend computation general options
%   The value of each property estimatedTrendE, estimatedTrendN and 
%   estimatedTrendV, as the trend estimation is carried out, is an 
%   'estimatedTrend' struct variable (see below). 
%   The value of estimatedTrendOpts is OptsGen.Trend, with OptsGen as
%   above.
%
%     detrendedZeroMeanE: detrended zero mean East time series
%     detrendedZeroMeanN: detrended zero mean North time series
%     detrendedZeroMeanV: detrended zero mean Vertical time series
%
%            CMEfiltered: CME-filtered time series and data
%                   CMEfiltered.t
%                   CMEfiltered.E
%                   CMEfiltered.N
%                   CMEfiltered.V
%                   CMEfiltered.estimatedTrendE
%                   CMEfiltered.estimatedTrendN
%                   CMEfiltered.estimatedTrendV
%                   CMEfiltered.Method
%
% The fields of these processing-related struct variables, also depending 
% on the kind of struct variable, i.e. 'offsets', 'outliers' or 
% 'estimatedTrend', as well as on the user's choices, are taken from
% the generated JSON files. They can be: 
%
%                 t1: initial time
%                 t2: final time
%                  N: actual mumber of days (gaps are excluded from count)
%     gap_percentage: percentage of gaps 
%                  k: number of estimated parameters (it is not kappa!)
%               ln_L: minimum value of log-likelihood ln(L)
%                AIC: Akaike Information Criterion (AIC=2*k+2*ln(L))
%                BIC: Bayesian Information Criterion (BIC=k*ln(N)+2*ln(L))
%             BIC_tp: another BIC value (BIC_tp=k*ln(N/(2π))+2 ln(L))
%              BIC_c: another BIC value, with some extra-penalties
%           ln_det_I: 
%
%         NoiseModel: struct variable depending on the user's choices:
%                     (if white noise is chosen): 
%                     NoiseModel.White 
%                                       NoiseModel.White.sigma
%                                       NoiseModel.White.fraction
%                     (if Powerlaw noise is chosen): 
%                     NoiseModel.Powerlaw 
%                                       NoiseModel.Powerlaw.sigma
%                                       NoiseModel.Powerlaw.d
%                                       NoiseModel.Powerlaw.kappa
%                                       NoiseModel.Powerlaw.fraction
%                     (if GGM (Powerlaw) noise is chosen):
%                     NoiseModel.GGM  
%                                       NoiseModel.GGM.sigma
%                                       NoiseModel.GGM.d
%                                       NoiseModel.GGM.kappa
%                                       NoiseModel.GGM.x1_phi  
%                                       NoiseModel.GGM.fraction
%                     (if FlickerGGM noise is chosen):
%                     NoiseModel.FlickerGGM  
%                                       NoiseModel.FlickerGGM.sigma
%                                       NoiseModel.FlickerGGM.d (always 0.5) 
%                                       NoiseModel.FlickerGGM.kappa (always -1)
%                                       NoiseModel.FlickerGGM.x1_phi  
%                                       NoiseModel.FlickerGGM.fraction
%                     (if ARMA is chosen):
%                     NoiseModel.ARMA
%                                       NoiseModel.ARMA.sigma
%                                       NoiseModel.ARMA.AR
%                                       NoiseModel.ARMA.MA
%                                       NoiseModel.ARMA.d
%                                       NoiseModel.ARMA.fraction
%                     (if ARFIMA is chosen):
%                     NoiseModel.ARFIMA
%                                       NoiseModel.ARFIMA.sigma
%                                       NoiseModel.ARFIMA.AR
%                                       NoiseModel.ARFIMA.MA
%                                       NoiseModel.ARFIMA.d
%                                       NoiseModel.ARFIMA.fraction
%
%         driving_noise: driving noise standard deviation (mm)
%                 trend: estimated trend (mm/y)
%           trend_sigma: estimated trend standard deviation (mm/y)
%
%                Sa_cos: yearly cos coefficient 
%          Sa_cos_sigma: yearly cos coefficient standard deviation
%                Sa_sin: yearly sin coefficient 
%          Sa_sin_sigma: yearly sin coefficent standard deviation
%          Sa_amplitude: yearly oscillation amplitude
%    Sa_amplitude_sigma: yearly oscillation amplitude standard deviation
%              Sa_phase: yearly oscillation phase
%        Sa_phase_sigma: yearly oscillation phase standard deviation
%
%               Ssa_cos: half-yearly cos coefficient
%         Ssa_cos_sigma: half-yearly cos coefficient standard deviation
%               Ssa_sin: half-yearly sin coefficient 
%         Ssa_sin_sigma: half-yearly sin coefficient standard deviation
%         Ssa_amplitude: half-yearly oscillation amplitude
%   Ssa_amplitude_sigma: half-yearly oscillation amplitude standard deviation
%             Ssa_phase: half-yearly oscillation phase
%       Ssa_phase_sigma: half-yearly oscillation phase standard deviation
%
%          jumps_epochs: []
%           jumps_sizes: []
%          jumps_sigmas: []
%
% See also tsDataFileIn, StaVelMain, geneOpts, InspOffsetMom, 
%   InspectOffsetSearch. 

% G. Teza, 2022

classdef tsData
    
    properties
                
        statName        % station name
        dateS           % date (string YYMMMDD) 
        dateyfrac       % date (fractional year)
        MJD             % modified julian date  
        GPSweek         % GPS week
        GPSday          % GPS day
        t               % date (MATLAB serial) 
        statLat         % station latitude (single value, degrees)
        statLon         % station longitude (single value, degrees)
        statHeight      % station heigh (single value, m)
        reflon          % reference meridian longitude (DEGREES)
        E0              % eastings (m), integer portion (from ref. meridian)
        Ed              % eastings (m), fractional portion
        E               % eastings (m), complete
        N0              % northings (m), integer portion (from equator)    
        Nd              % northings (m), fractional portion
        N               % northings (m), complete
        V0              % vertical (m), integer portion (from equator)    
        Vd              % vertical (m), fractional portion
        V               % vertical (m), complete
        antH            % antenna height (m) assumed from Rinex header   
        sE              % east sigma (m)
        sN              % north sigma (m)
        sV              % vertical sigma (m)
        cEN             % east-north correlation coefficient
        cEV             % east-vertical correlation coefficient
        cNV             % north-vertical correlation coefficient
        Lat             % latitude time series
        Lon             % longitude time series
        Height          % height time series 
        
        cleanedE        % cleaned time series - E component
        cleanedN        % cleaned time series - N component
        cleanedV        % cleaned time series - V component

        offsetsE        % offsets East
        offsetsN        % offsets North
        offsetsV        % offsets vertical
        offsetsOpts     % offsets general options
        
        outliersE       % outliers East
        outliersN       % outliers North 
        outliersV       % outliers Vertical
        outliersOpts    % outliers general options
        
        estimatedTrendE % estimated trend East 
        estimatedTrendN % estimated trend North 
        estimatedTrendV % estimated trend Vertical
        estimatedTrendOpts % estimated trend general options
        
        detrendedZeroMeanE % detrended zero mean E
        detrendedZeroMeanN % detrended zero mean N
        detrendedZeroMeanV % detrended zero mean V
        
        CMEfiltered     % for CME-filtered data
        
    end
    
    methods
        
        function obj = tsData(filena)
            
            if (nargin < 1) || isempty(filena)
                [fily1,pathy1] = uigetfile(...
                    {'*.tenv3';'*.tenv';'*.kenv';'*.pos'},...
                    'TIME SERIES ASCII env FILE (*.tenv3/*.tenv/*.kenv/*.pos)');
                if isnumeric(fily1) || isnumeric(pathy1)
                    filena1 = []; 
                else
                    filena1 = fullfile(pathy1,fily1);
                end
                filena = filena1; 
            end
            if ~isempty(filena)
                [~,~,fi] = fileparts(filena);
            else
                fi = [];
            end
            if contains(fi,'tenv3')
                opt = 4;
            elseif contains(fi,'tenv')
                opt = 3;
            elseif contains(fi,'kenv')
                opt = 2;
            elseif contains(fi,'pos')
                opt = 1;
            else
                opt = 0;
            end
            if ismember(opt,2:4)
                fid = fopen(filena);
                fgetl(fid);     % to remove header
                
                if opt == 4     % tenv3 ASCII file

                    C = textscan(fid,...
                        '%s %s %f %d %d %d %f %f %f %f %f %f %f %f %f %f %f %f %f %f');
                    Ccontrol = C{20};
                    if sum(isnan(Ccontrol)) > 2 % in this case, a new textscan is necessary
                        Iadd = true;
                        fclose(fid);
                        fid = fopen(filena);
                        fgetl(fid);
                        C = textscan(fid,...
                            '%s %s %f %d %d %d %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f');
                    else
                        Iadd = false;
                    end
                    C1 = C{1};
                    obj.statName = char(C1(1));
                    obj.dateS = char(C{2});
                    % C{3} not used because of problem of unadequate precision (possible double dates)   
                    obj.MJD = C{4};
                    obj.t = MJD2serial(C{4});
                    obj.dateyfrac = MJD2frac(C{4}); % to solve the problem of C{3} 
                    obj.GPSweek = C{5};       
                    obj.GPSday = C{6};
                    obj.reflon = C{7};
                    e0e = C{8};     % extracted value          
                    ede = C{9};
                    ee  = e0e+ede;
                    obj.E0 = e0e;          
                    obj.Ed = ede;          
                    obj.E  = ee;  
                    n0e = C{10};          
                    nde = C{11};
                    ne = n0e+nde;
                    obj.N0 = n0e;            
                    obj.Nd = nde;         
                    obj.N  = ne;
                    v0e = C{12};          
                    vde = C{13};
                    ve = v0e+vde;
                    obj.V0 = v0e;          
                    obj.Vd = vde;         
                    obj.V  = ve;           
                    obj.antH = C{14};         
                    obj.sE = C{15};          
                    obj.sN = C{16};          
                    obj.sV = C{17};          
                    obj.cEN = C{18};         
                    obj.cEV = C{19};         
                    obj.cNV = C{20};
                    if Iadd
                        Lat = C{21};
                        obj.Lat = Lat;
                        Lon = C{22};
                        Ineg = Lon < 0;
                        if sum(Ineg) > 0
                            Lon(Ineg) = 360+Lon(Ineg);
                        end
                        obj.Lon = Lon;
                        Height = C{23};
                        obj.Height = Height;
                        if numel(Lon) > 100
                            statLat = mean(Lat(end-100:end));
                            statLon = mean(Lon(end-100:end));
                            statHeight = mean(Height(end-100:end));
                        else
                            statLat = mean(Lat);
                            statLon = mean(Lon);
                            statHeight = mean(Height);
                        end
                        obj.statLat = statLat;
                        obj.statLon = statLon;
                        obj.statHeight = statHeight;
                    end
                
                elseif opt == 3     % tenv
                    
                    C = textscan(fid,...
                        '%s %s %f %d %d %d %f %f %f %f %f %f %f %f %f %f');
                    C1 = C{1};
                    obj.statName = char(C1(1));
                    obj.dateS = char(C{2});
                    obj.MJD = C{4}; 
                    obj.t = MJD2serial(C{4});
                    obj.dateyfrac = MJD2frac(C{4});         
                    obj.GPSweek = C{5};        
                    obj.GPSday = C{6}; 
                    obj.Ed = C{7};
                    obj.E  = C{7}; 
                    obj.Nd = C{8};  
                    obj.N  = C{8}; 
                    obj.Vd = C{9};
                    obj.V  = C{9};
                    obj.antH = C{10};         
                    obj.sE = C{11};          
                    obj.sN = C{12};          
                    obj.sV = C{13};          
                    obj.cEN = C{14};         
                    obj.cEV = C{15};         
                    obj.cNV = C{16};
                    Iadd = false;
                
                elseif opt == 2            % kenv
                    
                    C = textscan(fid,...
                        '%s %d %d %d %d %d %d %d %f %f %f %f %f %f %f %f %f');
                    C1 = C{1};
                    obj.statName = char(C1(1));
                    MJDe = C{3};
                    obj.MJD = MJDe;
                    obj.dateyfrac = MJD2frac(MJDe);
                    yd = C{4};
                    md = C{5};
                    dd = C{6};
                    dnum = datenum(yd,md,dd);
                    [w,d,~,~] = gpsweek(dnum);
                    obj.GPSweek = w;        
                    obj.GPSday = d;
                    obj.t = MJD2serial(MJDe);
                    obj.Ed = C{9};
                    obj.E  = C{9};
                    obj.Nd = C{10};
                    obj.N  = C{10};
                    obj.Vd = C{11};
                    obj.V  = C{11};
                    obj.sE = C{15};          
                    obj.sN = C{16};          
                    obj.sV = C{17};
                    Iadd = false;
                end
                fclose(fid);
                if ~Iadd
                    [Lat,Lon,ellH] = GetLonLat(obj.statName); % data from NGL
                    obj.statLat = Lat;
                    obj.statLon = Lon;
                    obj.statHeight = ellH;
                end

            elseif opt == 1     % case .pos data
                
                nr = linecount(filena);
                % header management:
                fid = fopen(filena);
                tline = fgetl(fid);
                nob = 0;
                nh = 1;
                while nob == 0
                    if numel(tline) > 8
                        tentdate = tline(1:8);
                        tentdaten = str2double(tentdate);
                        if isempty(tentdaten)
                            tentdaten = NaN;
                        end
                        if ~isnan(tentdaten)
                            nob = 1;
                        end
                    end
                    if numel(tline) >= 20
                        if strcmp(tline(1:14),'4-character ID')
                            statname = tline(17:20);
                        end
                    end
                    if nob == 0
                        tline = fgetl(fid);
                        nh = nh+1;
                    end
                end
                nh = nh-1;
                fclose(fid);
                % data extraction:
                fid = fopen(filena);
                tline = fgetl(fid);
                M = zeros(nr-nh,24);
                for k = 1:nr
                    m = k-nh; 
                    if k > nh
                        tline = erase(tline,'rapid');
                        tline = erase(tline,'final');
                        tline = erase(tline,'suppl/supplf');
                        tline = erase(tline,'campd');
                        tline = erase(tline,'repro');
                        tlinen = str2num(tline); %#ok<*ST2NM> 
                        M(m,:) = tlinen;
                    end
                    if k < nr
                        tline = fgetl(fid);
                    end
                end
                fclose(fid);
                %
                obj.statName = statname;
                mjd = round(M(:,3));
                obj.MJD = mjd;
                obj.dateyfrac = MJD2frac(mjd);
                td =  MJD2serial(mjd);
                obj.t = td;
                [W,D,~,~] = date2gpsw(td);
                obj.GPSweek = W;        
                obj.GPSday = D; 
                obj.Ed = M(:,17);   % note NEU sequence for .pos data
                obj.E  = M(:,17); 
                obj.Nd = M(:,16);  
                obj.N  = M(:,16); 
                obj.Vd = M(:,18);
                obj.V  = M(:,18);
                obj.sE = M(:,20);          
                obj.sN = M(:,19);          
                obj.sV = M(:,21);          
                obj.cEN = M(:,22);         
                obj.cEV = M(:,24);         
                obj.cNV = M(:,23);
                Lat = M(:,13);
                Lon = M(:,14);
                Height = M(:,15);
                obj.Lat = Lat;
                obj.Lon = Lon;
                obj.Height = Height;
                if numel(Lon) > 100
                    statLat = mean(Lat(end-100:end));
                    statLon = mean(Lon(end-100:end));
                    statHeight = mean(Height(end-100:end));
                else
                    statLat = mean(Lat);
                    statLon = mean(Lon);
                    statHeight = mean(Height);
                end
                obj.statLat = statLat;
                obj.statLon = statLon;
                obj.statHeight = statHeight;

            else
                warning('All properties are undefined');
            end
        end

    end

end